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Abstract. - We study acoustic wave propagation in a uniform stationary flow. We develop 
a method founded on the Lorentz transform and a hypothesis of irrotationality of the acoustic 
perturbation. After a transformation of the space-time and of the unknown fields, we derive 
a system of partial differential equations that eliminates the external flow and deals with the 
classical case of non advective acoustics. A sequel of the analysis is a new set of perfectly 
matched layers equations in the spirit of the work of Berenger and Collino. The numerical 
implementation of the previous ideas is presented with the finite differences method HaWAY on 
cartesian staggered grids. Relevant numerical tests are proposed. 
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1) Introduction 

The acoustic dimensioning of a civil aircraft requires the use of numerical models to predict 
the radiated acoustic field emitted by the engine. We study this problem in the case of 
advective acoustics. The specific case of classical acoustic wave propagation can be viewed 
as a scalar version of Maxwell equations. We use our previous experience acquired in the 
simulation of the propagation of electromagnetic waves using the finite differences method 
|DDS94j to rapidly develop a three-dimensional software simulating the propagation of 
acoustic waves. 

Our work is structured as follows. In order to take into account the external aerodynamic 
flow, we first come back to the equations of gas dynamics and consider the acoustic field 
as a first order linear perturbation of such a flow. Then we use physical ideas based on 
the Lorentz group invariance in Section 3, in the spirit of [AB86] . to deal in Section 4 with 
the case of advective acoustics in the same way as the non-advective ones. We develop 
in Section 5 a smart solution of the difficult problem of absorbing layers. The numerical 
aspect with the use of staggered grids is tackled in Section 6, and relevant physical and 
numerical tests are proposed in Section 7 

2) Non linear acoustics 

2-1 Barotrope gas dynamics 

We consider the propagation of sound waves in a uniform two-dimensional subsonic flow 
of a compressible fluid. This phenomenon is described by the nonlinear Euler equations 
for gas dynamics, see |LL54] for example, as : 



where u = [u, v) is the velocity vector, p the density of the fluid, p the pressure of the 
fluid and s the entropy. We also know that : 




< 





ds ^ ds , ds 



(2) 




where Cy is the calorific capacity at constant volume and {po,Po) a state of reference. 
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2-2 Linearization around a stationary state 

We linearize the system ([T]) around a constant state Wq defined by : 

Wo = {po,Uo,Vo,SoY , 

where (...)* is the transpose of a vector. The global state W of the system is defined 
around the state Wq thanks to the perturbation W = {p, u, v, s)* , as : 

W = Wo + W . 

A first idea of our approach is to use the impulses : 

pit = (po + p)(mo + u) = PoUq + ^ + pu 
pv = (po + p)ivo + v)= povo + C + pv , 

and to linearize them considering the variables p, u, v, s as first order infinitesimal quan- 
tities. We then introduce the linearized impulses : 

^ = Pou + puo 
( = PoV + pvo . 

We have the following classical hypothesis, see |LL54] : 

Hypothesis 1 Isentropy of the flow. 

The linearization of the fourth equation of the system ([T]) gives : 

ds ds ds _ ds ^ 
dt dx ^'^ dy dt 

If we consider the perturbation of entropy at the initial time to be null, that is to say 
s{x, y,t = 0) = 0, we deduce that s{x, y,t) = during the time evolution. 

Then the system ([1]) can be shared, first into a stationary aerodynamic system : 

div (poUo) = 
PoUo • Vuo + Vpo = , 

and then into an isentropic acoustic system : 

dt dx dy 

(3) { 1^ + — (2noe+ V2 V ) +— KC + ^o^-P^ot^o) = 



dt dx \ Cq / dy 

dC d d ( (? — v^ \ 

with p = CqP and Cq the speed of sound, deduced from the linearization of 
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Proposition 1 (Advection of the acoustic vorticity). 

If the external flow Wq = {po, Uq, Vq, Sq)* is stationary and uniform, the acoustic vorticity 

du dv 1 1 , n doj 

uj = — — IS advected by the flow, i.e. — = 0. 

oy ox at 

Proof. This property is classical, see |LL54] for example. We give the proof for complete- 
ness. We have from the system : 

d d dp 

= ^ + - PUo) + Vo^i^ - pUo) + Uo— + Uq— + — 

ot ox oy ox oy ox 

d d d dp 

dC d d f c^ — v"^ \ 

d( d ,^ , d ,^ , d( dp 

ot ox oy ox oy oy 

, d d d \ dp ^ 

We differentiate the flrst set of equations by y and the second by x, we eliminate the 
pressure field and obtain : 

1 d f d ,^ X d , ^ A d ,du dv , du 

^d^: " 'd-x^^' '"'^ ) = di^d-y-d-J = Tt='- 

□ 

Hypothesis 2 Irrotationality of the acoustic vorticity. 

If we consider the acoustic perturbation at the initial time to be irrotational, i.e. 
rot u (x, t = 0) = UJ (x, t = 0) = 0, we then deduce with equation (jl]) that rot u (x, t) 
during the time evolution. 



3) Lorentz Transform 

We consider the two-dimensional equations of advective acoustics when the velocity of 
the fluid is parallel to a particular direction; we suppose specifically : 

(5) u = Uo ex . 

We search a space-time transform (x, t) i — > {x',t') so that in the new space-time {x',t'), 
the pressure fleld is the solution of the wave equation. We flnd that this space-time 
transform is a Lorentz transform. With it, we derive a new set of equations and 
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prove that the corresponding system can be reduced to the classical case of non advective 
acoustics. 



3-1 Change of space-time 

Considering a flow of velocity given by equation the system is written as : 



(6) 



9i^d ( , , (eg - ul) 

at ox \ Cq 

dC ( C) 

. dt dx dy 



d_ 

dy 



(uoC) = 



which is a pleasant conservative form. We easily deduce that the pressure field p{x,y,t) 
is solution in the (initial) space-time {x, y, t) of a wave equation : 



(7) 



where A 



CqAp = 



+ Un 



+ — r IS the laplacian m two dimension space. 



dx'^ dy'^ 



Proposition 2 (Lorentz transform and equation of pressure). 

We suppose that the advective velocity satisfies equation (j5]). We define the Mach number 
as Mq = ^ and the Lorentz space-time transform as : 



(8) 



1 



X 



y' = y 
t' = t + 



X 



Mn 



Co(l - 



X . 



In this new space-time, the pressure field is considered as a function of the new set of 
space-time coordinates {x',y',t'), i.e. : 



(9) 



p {x ,y',t') =p{x,y,t) , 



and is the solution of the wave equation with a modified celerity : 



(10) 



d^p' 



,2 ^0 



dil - Ml 



- 



+ 



dx' dy' 



. 



reduced from the "pure" sound celerity by a similarity factor a/1 — Mq . 



6 



Frangois Dubois, Eric Duceau, Frederic Marechal and Isabelle Terrasse 



Proof. We first explain the way we derive the Lorentz transform ([8]) to remove the ad- 
vective contribution 2uq-q^ in equation ([7]). In the new space-time {x',y',t'), we want 
the pressure field to be solution of the wave equation. We search the new space-time 
coordinates {x',y',t') as : 



(11) 



X = ax 

y' = y 

t' = t + Px 



The transformed equation ([7j) takes the form : 



dxdt 

d d \ 2 a 

+ -'^^ 



p{x,y,t) 



((1 + uo^r - cl/3') + 2« (%(1 + no(3) - /3c^) ^ 



2/ 2 2^, 

-a (co - Uq) 



dx 



2 '^0 r^ ,2 



ay 



p{x ,y ,t) 



Then we obtain : 



(12) 



((1 + UoPr - clP') 1^ + 2a K(l + uof3) - Pel) ^ 



2/ 2 2^ 

-a (co - Wo; 



,2 H 



2 o /2 



oy 



p{x ,y ,t) = 



The conditions upon a and /3 to find the wave equation are clear from equation (I12p : on 
the first hand no further crossed partial derivation between space and time, that is : 



(13) 



2a {uoil + UoP) - (3cl) = , 



and on the other hand equality of the coefficients of double derivations in space to have 
a laplacian operator invariant by rotation : 



(14) a\cl - ul) = cl . 

The unique solution (a,/3) of the previous 2x2 linear system (fT5]) -( |T^ is : 



Co 

a = — ; 
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and with this set of coefficients, the space-time transform (fTT]) is exactly equal to the 

system ([8]). Moreover, we remark that the coefficient of in equation (fT2|) is now equal 
to : 

{{1 + uofiY - 



Cq Uq 



and, in our transformed space-time {x',y',t'), the pressure field defined by the condition 
([9]) is the solution of the wave equation f fTOj) . □ 



). We have the following 



3-2 Change of unknown functions 

Let us apply the Lorentz transform ([8]) to the acoustic system 
proposition : 

Proposition 3 (New unknown functions for advective acoustic). 

We assume that the hj^othesis 2 of irrotationality of the acoustic vorticity is satisfied and 
that the advective velocity field is defined by equation After applying the Lorentz 
transform ([8j) and the following change of pressure and impulse functions : 



(15) 



the advective acoustic system 



P 
I 



p + 



1 



jl - M2) ^' 



Ml 



(16) 



I can be written as : 



dti 
dC 



dx' 








Proof. We ffist use the hypothesis of irrotationality of the acoustic vorticity in the third 
equation of the system ([6]) and obtain : 

aC ^ <9(^ - pup) ^d^_ uodp 
dx dy dy Cq dy 

Secondly we introduce the Lorentz transform ([8]) into the system ([6]). We have the fol- 
lowing transform of partial derivations : 

1 d 



d 




dx 




d 




f 


dy' 




d 


dt 


~ dt' ■ 



M2 dx' 



+ 



d 



clil 



M2) dt' 



8 



Frangois Dubois, Eric Duceau, Frederic Marechal and Isabelle Terrasse 



We then substract the first equation of the system 
one and, using the following notations : 



multiplied by ^ from the second 



(17) 



p'{x',y',t') =p{x,y,t) 
^\x',y',t')=^{x,y,t) 
C'{x',y',t') = C{x,y,t) 



we find 



r dp' 

di' 



Uo 



M2 dx' 



'I 



M2) dt' 



Mi 



di' 



+ 



'dy' 

Co - ul dp' 



d^ 

'dil^ ^l- Ml dx' ' l-M^dt' ' d^l-M^dx' 



dc , ^, de , (eg - ui) dp' 

Co dy' 



9f + 



. 



We gather the terms associated with the same operator of derivation : 



d_ 
dt' 
d 



p' + 



MoCo 



1 - Ml 



+ Co 



5f '^t^m; 

dC 



dt' 



+ (1-M' 



^^'-<>d-x 



d 



dy' 



p + 



dx' 
d 
)x 

MoCo 



i' 



,dC' 



'1 - M2 



(1 



Mi 



p' + 



MqCq 

--Mi 
= , 



and we substitute into the previous system the new unknown functions (p , ^ , C) introduced 
in the system (fTSj) . Then the system of equations (fT6l) is satisfied. □ 



Remark 4. 

The major consequence of propositions [2] and |3] is the following (operational!) remark. 
The resolution of the advective acoustic system is absolutly identical to the one obtained 
without advective flow, but with a propagation celerity scaled by a factor a/1 — Mq. 



4) Lorentz transform for mult i- dimensional flows 

In the previous section, we dealt with the case of a velocity field described by the equation 
([5]). In |AGH99) . Abarbanel et al consider a multi-dimensional flow as a one-dimensional 
flow, after a correct rotation of the studied medium by an angle 6 = tan~^(^) and 
considering the new velocity to be u^ew = a/uq -|- Vq. We observe that in order to study 
numerically the influence of the flow for acoustic propagation near objects, such an idea 
imposes a remeshing of the geometry for each change of the advective flow. In our opinion, 
this process is not compatible with the use of finite differences and with operational 
industrial constraints. We propose in this section to generalize the Lorentz space-time 
transform to a multi-dimensional flow and to extend our approach with the help of space 
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affinities to a multi-dimensional flow under the same hypotheses as before. In the next 
two paragraphs, we present the generalization of the Lorentz transform respectively to the 
two and three-dimensional cases. Only the two-dimensional case is proven in the present 
document. The proof of the three-dimensional case can be found in |Ma2k] . 

4-1 The two-dimensional case 

We now consider a subsonic uniform flow described by a velocity vector : 



u 



{Uo,Vo) . 



With such an external flow, the linearized isentropic Euler equations for advective acous- 
tics are : 



(19) 



dt ~'~ dx 



2^ 



Ua] \ d UqVq 

-P + ^{uoC + VoC ^P) 



dy 



'dt ^ a^^"""^ ^ """^ ~ ^ ~ ' °^ ^ 



-P 





. 



-0 ^2/ \ "-o 
Proposition 4 (Simpliflcation of the acoustic system). 

We generalize the Lorentz space-time transform when the external flow verifles (fT8|) . We 
introduce a new set of space-time coordinates as : 



(20) 



X 



t' = t 



: X 



72- 



Uo 



4(1 - Ml 



I ^0 



the Mach number as Mq = — a coupling coefficient a between the two cartesian 
coordinates as : 

(21) a = - 



UqVq 



1-^Jl-^ 



and the new unknown functions p, ^ et C defined by : 

1 



(22) 



p = p + 

c = 



1-M2 



{uoC + VoC) 




c^l-M^ 



cl 1-M2 



cl [ cl 1-Mi 



1-^ 



or I 
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Under the hypothesis 2 of irrotationahty of the acoustic vorticity, the new form of the 
acoustic system ( IT9|) is the following : 



^ f + ^o|7(?-«C) + cg|7(C-<) = 



(23) 



9t 
dV 



+ 1 - 



dx' 



dy' 





. 



Proof. First, by substracting with a correct coefficient the first equation of the system 
f fTOj) from the two others, we find : 



dp 

d . uq d 

[i - -3P) + 



dt 
d_ 

di 



(C 



^^0 ^ 

3^) 



dx 

d uqVq d 
(moC —P) + 



,9 UqVq 

P I + — (^o4 —P) 



dx 



dy 



dy 



-p 





. 



Using the hypothesis of irrotationahty of the acoustic vorticity, we have : 



im — —p) 



dy 



d_ 

dy 
d 



voiPou + puo) —p 



{PoVqv) 



d 



d_ 

dy 



(PoVqu) 



= ^{voC-^p). 
dx Cq 

The previous calculation gives the new system : 



f dp 

di 



^'"'dx^'^'dy 







d , ^ uo , d 



--0 



aitt-^!') + «:;h« + "oC + (i 
-(C \ . 

L dt cg^' dy 



dx 

^P) + TT- H + ^OC + (1 



M'o)p] 
M'o)p] 





. 



In the new space-time defined by the change of space-time (120|) . we have : 



dx 

d_ 
dy 

d_ 

L di ' 



d 



1 «o dx' 
1 d 



+ 



+ 



Uo 



d 



cg(l - M2) dt' 



d 



^ _ 4 dy' cl{l - M2) dt' 



"^0 



dt' 
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With the notation introduced in the system (fTTI) . the transformed equations take the 
algrebraic form : 



d_ 

dt' 



Mi 



d 



1-^ 



d_ 
d_ 



1-^ 
1-M2 



72" 





2 9 



= 



1 _ ''o 
i 



72- 





d 



1-M2 



By the change of unknown functions 

d . 



(24) 



dt' 



+ c 



the previous system becomes : 







1 ^0 



dti 



dp 





. 



. 



We focus here on the fact that the pair (^', C) is present in the first equation of whereas 
the new unknown functions are C, and C. Nevertheless with the last two equations of the 
system (122|) . we have the following calculation : 



r = Wi-^(e-«c) 



where a is defined by equation (12T]) . We then find the final form (l23l) of the system of 
advective acoustics in the new space-time {x',y',t'). □ 



4-2 The three-dimensional case 

The generalization to three dimension space can be done without any major difficulty. 
We have the following proposition proven in |Ma2k] : 

Proposition 5. 

We assume that the velocity of the external flow is given by : 



(25) 



Uo = {Uo,Vo,Wo) , 
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and that the Hypothesis 2 of irrotationahty of the acoustic vorticity is verified. We define 
the Mach number as Mq = ^ — - — - and a change of space-time 



as 



(26) 



X — 



z — 













1)2 
^0 







X 



t' = t + 



X + 



y + 



Wo 



c§(l - Mi) c§(l - Mi) ' cl{l - Mi) 



We introduce three couphng coefficients as : 




and a change of unknown functions as : 





«0 / 











1-M2 d I -Mi cl I -Mi 



C= 1/1 - T2 ^+ (1 12 )l IJ2 + 



Mi ' cl 'I- Mi cl I- Mi 

vqWo C ,^ ul + v^. x' 



2 d 1-Mi ' cl 'l-Mi 



The acoustic system takes the new form 
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dp 
di' 



< 



dti 
dt' 




ix) + c; 



dz' 



(x-/3e-7C) = o 



5) Acoustic Absorbing Layers 

Physical wave phenomena modelhng takes often place in the infinite two or three di- 
mensional space. Due to finite computing resources, the numerical simulations of such 
phenomena must be truncated to confined domains, then numerical artifical boundaries 
must be considered. Generally, numerical refiections of outgoing waves from the bound- 
aries of the numerical domain reenter the computational domain and falsify the results. 
Various methods have been proposed to reduce the infiuence of reentering waves in the 
computational domain. 

For many years, the numerical physicists, Israeli-Orsag |I081| . have developed the idea of 
layers of absorbing materials. Then the mathematical study of non refiecting boundary 
conditions has been developed after the pionnering work of Engquist-Maida |EM77j . The 
discrete studies of such absorbing conditions have been realized for scalar waves by Bayliss- 
Turkel |BT80) . for electromagnetic waves by Joly-Mercier |JM89j . Tafiove |Ta98] and for 
sismic waves by Halpern-Trefethen |HT86| among others. 

The current perfectly matched layers approach has been introduced by Berenger |Be94] 
in the context of computational electromagnetics; a mathematical interpretation of this 
model has been made by CoUino [Co85j . In |Hu96] . Hu proposes an adaptation of 
Berenger's model for advective acoustics. Nevertheless, Abarbanel et al |AGH99j and 
Rahmouni |Rah01| and |Rah99j have demonstrated that this model is mathematically ill- 
posed, i.e. that, if truncated to the first order terms, there exists a perturbation as small 
as we wish that can make the model unstable. These authors also propose well-posed 
models. 

Our approach uses a model of the type introduced by Hu. We focus in our study on 
the practical disadvantages to deal with a mathematical model whose principal symbol 
corresponds to a ill-posed problem. 
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5-1 Acoustic absorbing layers without external flow 

We propose in the following a precise description of acoustic absorbing layers, and we 
follow the ideas developed by CoUino |Co85] . We consider a semi-infinite medium in the 
x-direction defined hj Q = U Q~, where : 

(^Y) ( n- = y<0} 

^ ' I n+ = {ix,y), 0<y<5}, 

Q'^ representing the absorbing layers domain. 




Absorbing 
layers 



y = 5 



y = 



Physical 
domain 




Proposition 6 (System of acoustic absorbing layers). 

A system of partial differential equations that models absorbing layers of acoustic waves 
in the domain introduced in (127|) can be given as : 



(28) 



'Py 







where the absorbing coefficient cr*(y) satisfies : 

(29) [0, 5]3y\ — > a*{y) G R+ , a*(2/) > if y > and a*(0) = . 

The pressure p is defined hy p = p^ + Py . 



Proof. We follow essentially the idea of Collino |Co85| . The main idea to establish the 
acoustic system inside the absorbing layers is to introduce the Fourier-Laplace transform 
and to write the system fll9p in the complex plan. 
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We define the Fourier-Laplace transform by : 

v{h,y,u) = j jv{x,y,t) e-^(^*+'=-^) dxdt . 
Then in the domain = {{x,y), y < 0}, the system (fT9|l takes the form : 



(30) 



tup = -ik^CQ^ ~ 



iuj^ = —ikxP 
dy 



and the solution of the system (130|) for a propagation in the growing y-direction is obtained 
after the integration of an ordinary differential equation of degree 2 : 



(31) 



f 

V 


= poe-^'yy 




< i 


_ , ikyy 


with 


CO 




c 

v 


= -P^kye-^'yy 

CO 





~ki 



We establish a modified form of the system (130|) in that ensures that waves leaving the 
domain are not reflected back. Let us extend the variable y in the complex plan by adding 
an imaginary part depending on the function a* defined by equation (12^ . and that equals 
zero for a* = 0. We write precisely in the complex variable y parameterized by a real 
variable r], < r] < 6, as : 



y = v^iv) = V + 



o*{u] dw . 



We can draw the function y = (f{f]) in the complex plan as : 

Re y 

5 




For V equal to one of the variables of pressure or momentum, we introduce the function 
v{y), with y = ^9(77), as a function v* of the real variable r] : 



v*{v) = , V e{p,^,C} 

An elementary calculation gives us : 
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v*{r]) 



V ( v -\ / (T*(u) du] 

\ Jo J 



Vq exp 



7] + 



a*(u) du 



V (jj) exp ^ — - j a*{u)du 



We then deduce the following important property : 



(32) 



\v*{r,)\ < \viri)\ , 0<r]<5 , ve{p,^,C} ■ 



The property (l32l) is a consequence of an exponential decay of all the variables inside 
the absorbing layers. It is possible to derive the system of partial differential equations 
satisfied by those fields. A first algebraic calculation gives us : 



di) dv* drj 



dv* 



dy dr] dy iu + a* {rj) drj 

then we obtain with this new set of unknown functions a system in the {x, rj) domain 
issued from (150]) : 



(33) 



2t* 



100 



dC 



iuj + a*{rj) drj 



lUJ 



icu + a*{rj) drj 

The first equation can be rewritten while splitting the pressure field into two sub-pressure 
fields as : 

P* =P*x+pI , 

with p* and p* solutions of : 



lup^ 



,2f* 



iuj 



,dC 



^ iu + a*{rj) ° drj 
Taking the inverse Fourier-Laplace transform of the new system, we obtain ([28 



□ 

Considering a square domain [0, L]^, we define the thickness of the absorbing layers 
by 5x in the x-direction and 5y in the y-direction. The interesting studying medium is 
then [5xiL — 5x] x [5y, L — 5y]. We have the following proposition that generalizes the 
proposition [6] : 
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Proposition 7 (General acoustic system to solve). 

We consider two smoothing functions a* and a* defined by : 

a*{x) > Oifx e [0 ,6x[x] L - 5x, L 



(34) [0, L]3 XI — > a*(x) G R+ 

(35) [0,L]33y^a;{y)eR, 



cr*(x) = ifx e [S^,L-6cc], 

a;{y)>0ifye[0,dy[x]L-6y, L] 
a;{y)=Oiiye [6y,L-6y]. 



The acoustic system in the studying medium and in the absorbing layers can be written 
as : 



(36) 



dp. 

^ , , , d , . 

I dt + ""y^y^^ + d^(P^ + Py) 











Remark 9. 

This set of equations is the same as obtained by Hu in |Hu96] using velocity fields rather 
than impulses. 



Proof. The proof is similar to the one of proposition |6] 



□ 



5-2 Plane wave analysis 




The acoustic system, inside the absorbing layers in the y-direction is, after a Fourier 
transform : 

iujpx - iKcl^ = 

iujpy + a*{y)py-ikyclC = 

iuC, — ik,j.p = 

iu( + (J*{y)C — ikyp = . 
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The incident wave is a solution of the system : 



(37) 



iup - iklcli 



iuj^ — ik^p = 
iu 

iuC - -. — ■ — -iklp = 



iu + al ^ 

The reflected wave is a solution of the system : 

tUPr — ik^CQt,r — -. 

iu^r — iklpr = 



(38) 



-iklpr = . 



The transmitted wave is a solution of the system : 



(39) 



iupt - 




< iu^t - 


iklpt = 


iujQt - 


iklpt = . 



At y = 0, we write the continuity of the pressure field. We have : 
Pr = Rp and pt = Tp with 1 + R = T. We then have : 



= i?e and 6 = re 

(r = RC and (t = TC 



We deduce from the system ([3 



—p and C = -■ — 7—; 



p 



and we know that k^ = k], and ky 



-ky. We then have : 



C + Cr 



hi hi 
—p + —Pr 

uj . 00 



ik, 



UJ 



[I + R)p 



y 



iu + a 



7(1 



The system (139|) gives us : 



^, = ^{l + R)p = ^Tp 
u u 
ik:, ik:, 
6 = -^(1 + R)P = ^-^Tp 



lu + a\ 



iu + o\ 



P + Pr = Pt then 
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At the interface y = 0, we write the continuity of C. We then have Q 
write as : 



C + Crj that we 



(40) 



ikl 



;i + R)p 



ikl 



iu + alio) 



(1 - R)p 



5-3 Mathematical property of the absorbing layers system 

The system of acoustic absorbing layers ( 13 6 p can be written as : 



(41) 



dw ^dw aw ^ 

+ A—- + B—— + CW = {) 



dt 



dx 



dy 



where W = {p^, Py, ^, C)* 



A 



/ cl \ 



110 





B 



/ \ 

eg 



110 



and 



C 



o;{y) 










\ 




The principal symbol, M = ik^A + ikyB, of the system (HTj) is given by : 







/ ik^cl \ 



(42) 



M 



ikycl 






\^ iky iky 

The eigenvalues and eigenvectors of the principal symbol are re-evaluated in table 1. 



Eigenvalue 



Eigenvector 



(double) 




(- 


1,1,0,0)* 




CQsJ-kl - kl 


ikr^c^^ k^. 




"^k^^^o^J —kl — k^ 


h 

\ y Y 

kx 


^ ^2 _|_ ^2 




kx{kl + kl) ' 


-co^-kl - kl 


AkxCQ^J k'^ 


Ky 


^kyCo^y—k'^ — k^ ^ 


h 

'^y y 


^ ^,2 1 Z,2 

/Vj, Hy 




kx{k1 + ky) 


kx 



Table 1. Eigenvalues and eigenvectors of matrix 
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We notice that is an eigenvalue of muhiphcity order equal to 2 associated with a one- 
dimensional eigensubspace. The system (141 p is not hyperbolic and classical results relative 
to well-posedness of such systems (see |Rau91] ) can not be applied to the absorbing layers. 
There exists an arbitrarily small perturbation of the Cauchy problem for the system fl36l) 
with a*{x) = cr*{y) = that makes the system (HT]) ill-posed for or Sobolev norms of 
order 1. Nevertheless, our choice of the system 0411) does not produce unstable numerical 
results. 



Proposition 8. 

If we look for a solution of the form W 



Lp{t) e ^''^^e ^''^'^ of the system (HT]) . where 



Vjvf2 = (0, 0, ky, —k^y, the scalar function Lp{t) is an exponential decay in time. 

Proof. To establish this result, we determine the characteristic subspace for the eigenvalue 
A = 0, i.e. we calculate ker{M'^). We have : 



(43) - 



k 



2 r, 2 

y '^y 




V 















kx 

kxky 






kx ky 
ky^ 



\ 



and ker(M2 



/ 1 \ 

-1 


LV / 



/ \ 



V 



/J 



We analyse the stability of the system (HT]) under a perturbation following the direction 
of the eigenvector of that is not eigenvector of M. We note Vj[p = (0, 0, ky, —k^Y 
this vector that is a simple impulse pertubation in the direction orthogonal to the wave 
vector. We choose a state vector W of the form : 



W = ip{t) e-'^^'^e-'^yy Vm2 . 



In this case, the system (l36l) is written as : 



(44) 
with 



MVm2 



\ 



Vy 

~7ikrf. k' 





y 



CVm^ 



J 



\ 



\ 



a*{x)ky 
-(^*y{y)kx ) 



We then deduce that the first two equations of (jH]) impose that k^ky 
kx = and ky ^ 0, the third equation of (jH]) gives us : 



0. Therefore, if 



(45) 



'dt 



ky + al{x)kyip = i.e. 



-Q^ + <^x{^)v = ^ 
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We assume that cr*(x) > inside the absorbing layers, the solution of the equation (H5|) 
is an exponential decay in time of the function ip. □ 



We obtain the same result considering ky = and k^. ^ 0. Even if the principal symbol of 
the system (1411) is associated to a "ill-posed mathematical problem", the form of the zero 
order terms shows that even exciting the absorbing layers system in the direction of the 
characteristic vector, the perturbation is dissipated. 

We would like to predict the behavior of our absorbing layers model. Thus, we simplify 
our set of equations to the simplest model and study it. We establish the following 
proposition : 

Proposition 9. 

We consider the simplest 1-D non-hyperbolic model inside the absorbing layers, excitated 
with a source term centered at {xa,ya), in the direction of the eigenvector Vm^- We 
note W = {u,vy the state vector, Sx^^y^ the Dirac mass at the position {xa,ya) and we 
assume that the coefficients ai and a2 are strictly positive. The problem : 



(46) 



dW 
~dt 



+ 



1 




dW 
dx 



+ 













0-2 



w 







is stable as long as -(/^(t) is bounded. 
Proof. The system (146|) is written as : 



(47) 



' du dv 

qt ox 

dv , / \ r 

— +a2V = ^{t) 5xa,ya 

^ u{0) = v{0) = . 



The solution in f is f 

^xa,yai ^6 then deduce that u is solution 

of: 



du 
'di 



V.(5)e--(-^)d^U;^,,„. 



Then, the solution in u is u{t) = ii(t) 6'^^ y^, where \imt_^oo l^{t) = 0. Then, for a function 
'il){t) bounded, m — )■ and v is bounded. The system (H7|) is stable as long as ai > and 

0-2 > 0. □ 
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This simple model makes us think that our absorbing layers model is stable even if we 
are in the worst situation, i.e. there is a source in the direction of the eigenvector Vjvfa. 
Numerical simulations proposed in Section 7 confirm this result. 

One challenge in the future is to theorically understand the behavior of the solution of 
the true problem in the absorbing layers defined by the following system : 



dW .dW dW 



(4J 



+ A^ + + CW = m 5^.,y. Vm^ 



dt ' dx 
t^(0) = 



Nevertheless, the qualitative behavior proposed for the system (146|) gives a good idea of 
the behavior of the system (1481) . see Section 7. 



5-4 Acoustic absorbing layers with external subsonic flow 
Proposition 10 (General absorbing layers in two dimensions). 

We assume that the velocity for the external subsonic flow is u = ('Uo,fo) • A system of 
partial differential equations that models absorbing layers of acoustic waves is given by : 



( dPx 



dt 



(49) < 



CqVo 



8t +c,^l-MSa-(y)p„+ ^^—^ 



dx 







cr*{y)C + c',-^ = 



dy 



d 



dp d 



^ + ^(2«oe + ^oC) + (l-Mo^)£ + -(«oC) + 



co(i + 



dx dy 



+ ^^^^^ ia%x)p. + a%y)p,) + -"-"^^^j^^ + --:(^)^C = 
Co 



dt 



+ ^^^^ (a-(x)p.. + a-(.)p.) + - , 

Co Co a/1 - 

where a*{x) et (J*{y) are smoothing functions defined by (135|) . 



Proof. Here are the main ideas of the proof; the details of all the calculus can be found in 
|Ma2k| . If we consider a two-dimensional flow, we have shown that the acoustic system in 
the new space-time defined by (120]) is given by (1231) after the change of unknown functions 
( I22p . Using the method described in the section, we easily show that a general system 
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of dimensionless partial differential equations for the acoustic absorbing layers can be 
written in (x', y', t') using the functions a*{x) and (J*{y) introduced by equation (135|) as : 



(50) 



dt 



I + co^l-Mia*{y)py + c^ — (C-a^) 












dC I ~ d 

y^, + c,^^-Mia\y)C + -^,{p. + Py) = ^ 



where a is a coupling coefficient given by : 



a 



UqVo 



4 V 4 



We now wish to write the system (I50p in the initial space-time (x, y, t), using the initial 
unknown functions p, ^, C- We have : 



(51) e - aC 

and we easily calculate : 



1 



■^0 



1)2 
"^0 



(52) 



d_ 
dt' 
d_ 

dx' 

d_ 
dy' 



d_ 
di 



1 - 



cl dx cg(l - Ml) dt 



2 Q ^OA/l-i d 



cl dy cl{l - M2) dt 
We substitute ([5l]) and ([52]) in the system ([50D : 



dt 



^ + co\l-Mia*{x)p, + cl\l 



1^(e-«c) 

Cn dx 



^^(e-«c) = o 



1-M2 dt 



dt 



y +c,Jl-Mia\y)py + cl.l 



Vox I - ^ Q ^ 



1-M2 at 



7^ 



at 



at 



. 
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Using the change of variables ( 122|) . we deduce : 

1 



(53) 



p = p + 




cl V cl l-Mo^ 



We substitute the change of variables (153|) into the previous system, we then obtain : 



(54) 



dt 



dt 



dt 



uq dp 



Cn Ot 



dp d 



Mi 



+ 



d( vq dp 
dt ~ 4di 



+ (i -0^ + ^Ke + ^oC) = o 

. coyrr^.*(.)(^^ + (i 



d 



dp 



+ 7rKe + ^oC) + (i-Mo^)-^ = o 



ay 



dy 



Mi 



Un 



cl'l-Mi 



We use as new unknowns : 



Wo 



we then have : 



and we finally obtain : 



M'^ 



P = Px+Py , P = Px+Py 



^ + Co\/l-Mia*{x)p^ + cl^ 



' dp . 



^ ' dx ' dx 

dC vodp r — . f uovo 



2^ 



) ^odp n^2^9p 



+ |Ke + ^^oC) + (l-Mo^)| = 0. 







r2 ^ 1 



c 
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We wish to find a dynamic system, then we ehminate the terms in the last two equations 
of the system (155|) . Using the equality p = Px + Py, and adding the first two equations of 
the system fl55p . we deduce : 

(56) ^ + + co^ + coVl - {y)p, + ^7^=^^ (^)C = , 

and we substitute the equation (15^ in the last two equations of the system (I55p . Hence, 
we find the result fl49l) that ends the proof. □ 



Remark 13. 

Various authors propose a system of partial differential equation for absorbing layers for 
advective acoustic (see |AGH99] . |Hu96] . |Rah99] . [RahOl] for example). Each of them 
have to solve 6 equations, whereas we propose a system composed by only 4 equations. We 
see this property as a consequence of our precise physical analysis based on the Lorentz 
transform and our change of unknown functions. 



6) Discretization with the "HaWAY" method 

This section deals with the numerical resolution of the equation of acoustic without an ex- 
ternal flow. We use finite differences with staggered grids as introduced by Harlow- Welsch 
(MAC method) |HW65j . Arakawa (C grids) [Ar66j and Yee |Yee66| for electromagnetism. 
HaWAY comes from Harlow- Welsch, Arakawa, Yee. The acoustic system can be written 
as : 



(57) 



dt' 


^""'dx 
dp' 


It' 


dx' 
dp' 


dt' 


dy' 



2dC 
^'dy' 







^ + =0 



-^ + ^ = 



We first propose to nondimensionalize the previous system. We then explain the numerical 
scheme chosen in the free space and in the acoustic absorbing layers. 



6-1 Dimensionlessness of the acoustic system 

This section is introduced for the completeness of our meaning. We refer to |Se75] for this 
kind of purpose. We nondimensionalize the set of equations (157|) by writing each variable 
as : X' = X*X, for X' pressure, impulse, time and space variables and X* a reference 
dimension for each variable: a reference pressure p*, reference impulses ^* and (*, a time 
reference t* and reference lengths x* and y*. The new form of the system (15 7p is : 
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t* C 2 C 9C 



~'~ ^ p* X* dx ~'~ ^ p* y* dy 



dp 
di 

^ ^i^dp ^ ^ 

dt ^* X* dx 
dC ^ r_p^dp_ ^ ^ 

, dt C* y* dy 

We decide here to choose the followinEf coefficients Cn , Cn , and — — equal 



to 1. We then deduce that we have : 

1 



-p 



1 



X 

Co" " Co" t* 

and the resulting set of dimensionless equations is 



Co 



t* p* 



l£_ 

t* 



Co 5 



(58) 



dt dx dy 
df dp 



^ dt ~^ dy 



6-2 Staggered grids for acoustics 

By analogy with electromagnetism (see |DDS94| ). we decide to use the cartesian staggered 
finite differences method to solve the system 0581) . We decompose a model domain Q = 
[0, L]^ into finite elements with an isotropic meshing of space step Ax = Ay = j {J E N* 
is the number of cells in each direction). The cell K,,i , , i is defined as : 

=] lAx , {i + 1) Ax [ X ] J Ay , (j + 1) Ay [ 

We share the time with the help of a time step At and introduce the n*^ "entire time" 
t^ = nAt . By convention, we know that the pressure is defined at entire times t" in the 
center of the mesh K^_^_l j^i and that the impulses are defined at semi-entire times t"''^^ 
on the edge of the mesh. The variables in a mesh are defined as below : 
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The numerical scheme used in the free space is in two-dimensional space : 
• Discretization of the pressure equation : 



rP+i — rP 



C3) + (Ci 



CD 



• Discretization of the impulse equations : 



.n+f 
3 



where we have cr = 



At _ At 
Ax Ay 



This numerical scheme is an explicit second order in time and space scheme, stable under 
the Courant-Priedrichs-Lewy condition : 



(59) 



CFL I 



At < 
At < 



1 + 1 + 1 
Ax^ Ay^ Az^ 



in two dimension space, 
in three dimension space. 



The boundary condition is supposed to be on the edge of the mesh. The Dirichlet bound- 
ary condition is written as : u.n = . 

6-3 Numerical acoustic absorbing layers 

We deduce from section 4 the set of equations to solve in the absorbing layers : 







— + a {x)p^ + — 



Q^+<y*{y)Q+Q^{px + Py) = Q, 

where (t*{x) and cr*(y) are the smoothing functions strictly positive inside the absorbing 
layers. The set of equations is ended by a Dirichlet condition on the edge of the whole 
studied domain : u.n = , where n is the external normal to the domain. The 
discretization of such a boundary condition for the whole studied domain is : 
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,0 



0, 0<i<J-l, n>0 



C4!j = 0' 0<^<J-1, n>0 



We propose to use the same discretisation as before. In the absorbing layers, we have to 
know ^ and ^ respectively at times n + letn + 1. We decide to center those 

values in time, we write : 



n+7 



<rn+i 



n+1 



2 V'^».i+I "^^'i+l 



2 v''»+ij 



/■■- ' 2 



The numerical scheme, while noting a — 



At _ At 
Ax ~ Ay 



(isotrop meshing) , can be written as 



[i + \) At 



2^ C 



2 + a*(« + i) At 
2-a;(j + i) At 



2 + fT*(j + 



2 
2 



1 

At 



At 



2 + a*(« + i) At 



2(7 



2 + ^;(j + i) At 



^ j) At +1 



2 + (7;(j) At 



2 + a*(«) At 

2(7 

2 + (7*(j) At 



with p = p^+py . 
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7) Numerical tests 

The computational domain is shared into two areas : the studying medium and the 
absorbing layers. Those two areas are defined as shown below : 



Xr 




Absorbing layers 



Studying medium 



7-1 Mathematical experiments 

We first propose to numerically validate our mathematical analysis of the absorbing layers 
acoustic system (HTj) with a test case proposed by O. Pironneau |Pi99] . We decide to place 
an acoustic pulse inside the absorbing layers, in the direction of the eigenvector associated 
to the eigenvalue in order to enforce the unstability due to the lack of hyperbolicity as 
shown before. For this test, the computational domain, symetric in the x-direction and 
the ^/-direction, is defined by Xmax = 5 and Lpmi = 45. We consider two different tests. 
For the first one, we solve : 



(60) 



+ A—- + B—- + CW 



dt 



dx 



dy 



For the second one, we solve : 



(61) 



+ A—- + B—— + CW = 



dt 



W{0) 



dx 
( 



dy 








exp(-(ln2) g 



exp f— (ln2) 



jx-Xg)^ +{y-ya)'^ 
9 



{x-Xa)^+{y-yaY 
9 



For the first test case, the excitation ip{f) Vm^ is centered inside the absorbing layers at 
i^a, Va) = (25, 0). The absorbing coefficients in the absorbing layers are constant in the x 
and in the y-direction. The reference excitation is : 
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ip^y = exp -(In 2) 



{X - Xaf + {y- Vaf 
- c^P \ -V^^^l ^ 

We solve the problem (l60l) with the particular data given by : 



PiXa,ya,'t) = 
^{Xaiyai^) 
CiXa,ya,t) = 



dx 

The observing points are centered at = (45,0), (2:2,2/2) = (25,0), (x^jy^) = (0,0), 

(x4,y4) = (—45,0), (2:5,2/5) = (0,25) and (xgjj/e) = (0,-25). We observe the results at 
{x2, 2/2) = (25, 0). We obtain the graph presented of Figure 1. 

-3 



X 10 



Pressure 

x-lmpulse 

y-lmpulse 




500 1000 1500 2000 2500 3000 

Time 

Figure 1. Pressure and impulse fields at {x,y) = (25,0). 



We remark that, for long times, the pressure field converges to and the impulse fields to 
a non-zero value. We obtain the results predicted by our simple 1-D model (see proposi- 
tion [9]) noting u the pressure field and v the impulse fields ^ or (. We notice that the sign 
is changing if we consider a positive or a negative source. We observe approximately the 
same results for the other observing points. The pressure field converges to except at 
(x5, 2/5) and {xq, 2/5), where there is a slight residual at (xs, 2/5) and —Vp at {xq, y^). The 
^ impulse field converges to at (xi, 2/1), [x^, 2/3) and {x^, 2/4). The residual at (X5, 2/5) is 
and — at (xg, yo). The ( impulse field always converges to a constant value as predicted 
by the 1-D model. 
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For the second test, we solve the problem (l6Ti) with an excitation centered inside the 
absorbing layers at {xa, Ua) = (25, 0). As before, we analyse the results at (x2, 2/2) = (25, 0). 
The evolution of the pressure and impulse fields are drawn in the following figure : 



0.02 
0.01 



0.01 
0.02 
0.03 
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y-lmpulse 
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3000 
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Figure 2. Pressure and impulse fields at (x, ?/) = (25,0). 



All the fields converge to for all the observing points. We obtain the same results as pre- 
dicted by the proposition [9] with an excitation function Tp{t) such as '0(t)dt is bounded. 



The conclusion of these mathematical experiments is that even if we excite this non 
hyperbolic system in the direction of the non caracteristic vector, the zero order damping 
terms insure that there is no numerical explosion of our results. 

7-2 Physical experiments 

We have proven the stability of our absorbing layers. We now want to study numerical 
reflections of outgoing waves from the boundaries of the computational domain for various 
speeds of the external flow. We then consider the problem (jH]) in the computational domain 
defined by x^ax = 25 with the absorbing layers outside {Lprni is now a parameter), an 
acoustic pulse centered at [xa, Ua) = (0, 0) and an excitation in the right hand side [p, ^, (Y 
given by : 

{x - Xaf + {y - VaY^ 



(62) 



p{x,y,t) = exp ( -(ln2)- 

^ix,y,t) = 
ax,y,t) = 0. 



9 



sin(7rt) 



We compare the calculated solution, denoted by p, to the numerical solution obtained in 
the domain defined by x^ax = 150. We take Ax = 1 and At following the CFL ( !59|) . For 
t < 300 At, it is easy to see that no reflection from the boundaries can interact with the 
solution within the small domain [—25, 25]^ and such a solution is the numerical solution 
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in an infinite domain for t < 300 At. This solution is considered as a reference, noted 
Pref, and the computation of the error \p — Pref \ for each time step indicates the efficiency 
of the absorbing layers. As explained in Section 5, the boundary condition is imposed on 
the edge of a cell and the pressure is calculated in the middle of the cell. The observing 
point is then taken only half a cell near the absorbing layers at (25,0). 




Figure 3. L2— error of the pressure for 4, 10 and 20 absorbing layers, ^ = (0.5, 0). 




50 100 150 200 250 300 

time 

Figure 4. L2— error of the pressure for 4, 10 and 20 absorbing layers. 



Co 



(0.5,0). 
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Figure 7. Iso-lines of the pressure field for t = 40 At, t = 80 At and t = 120 At. 
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^ure 8. Iso-lines of the pressure field for t = 40 At, t = 80 At and t = 120 At. 
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We first consider a domain without external flow. We notice that the L2-error of the 
pressure field p computed for each time step at the observing point (25, 0) for various 
thickness for the absorbing layers. We notice that for small absorbing layers (4 cells), we 
have a good accuracy for the results, and that increasing the thickness from 4 cells to 20 
cells improves the accuracy by 2 orders of magnitude. 

We also compare, for various thickness of the absorbing layers, the "exact" and the numer- 
ical solution considering the same observing point and the same acoustic source, for three 
velocity vectors defined by u = (0.5co,0), u = (^Cc^Cq) and u = (-^cc^^Cq). 
We notice an improvement of the accuracy by 2 orders of magnitude for an absorbing 
layers growing from 4 cells to a 20 cells. 

The results are satisfying. Nevertheless, when the number of cells in the absorbing layers 
is increasing, the error is small but remains measurable, even for the long times. We think 
that this behavior could be improved in future work. 

Conclusion 

We have explored a new method for solving the equations of advective acoustics based on 
a change a space-time variables (Lorentz transform) and a change of unknown variables. 
We have also derived a system of equations (H^ to modelize the absorbing layers for the 
acoustic model. The system of partial differential equations established in the absorbing 
layers is well-posed due to the zero order term. The staggered grid "HaWAY" method has 
been used for the numerical implementation and experiments have proven the efficiency of 
such a method. When we force a punctual acoustic source inside this numerical domain, 
our experiments show that the results remain bounded and our method is stable from 
a practical point of view. Notice that we explain our method only in two-dimensional 
space, but we extend it easily to three-dimensional space. 
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